#include "Tria3ThermalCalculator.h"
#include "Utility/Algebra/Tools.h"
#include "DBManager/IService/Input/Node/INodeService.h"
#include "DBManager/IService/Input/Element/IElementService.h"
#include "DBManager/IService/Input/Subcase/ISubcaseService.h"
#include "DBManager/IService/Input/Property/IPSHELLService.h"
#include "DBManager/IService/Input/Material/IThermalMaterialService.h"

#include "SGFem/Pre/Assembly/Assembly.h"
#include "SGFem/Pre/Assembly/ModelInfoQuery.h"


using SG::Algebra::Point;
using SG::Algebra::Vector3D;
using SG::Algebra::Matrix;
using SG::Algebra::identity;
using SG::Algebra::det3X3;
using SG::Algebra::length;
using SG::Algebra::calcTriArea;
using SG::DataStructure::Common::Real;
// using SG::DataStructure::Common::TOLERANCE;
using SG::DataStructure::FEM::PSHELLData;
using SG::DataStructure::Common::Id;
using SG::DataStructure::FEM::ThermalMaterial;
using namespace SG::DataStructure::FEM;
using namespace SG::DBManager;
using SG::FEMSolver::AssemblyEigen;
using SG::FEMSolver::Element::Tria3ThermalCalculator;
std::size_t Tria3ThermalCalculator::s_nNode = 3;
std::size_t Tria3ThermalCalculator::s_nDim = 2;


void Tria3ThermalCalculator::Initialize (SG::DataStructure::FEM::ElementTypes type)
{
    m_type = type; 
}

void Tria3ThermalCalculator::ComputeStrainMatrix (const Matrix& vertices, const Real &area,  _OUT Matrix& B) const
{
    //Real ai, aj, ak;
    Real bi, bj, bk;
    Real ci, cj, ck;

    //ai = vertices(1,0)* vertices(2,1) - vertices(2,0)* vertices(1,1);
    //aj = vertices(2,0)* vertices(0,1) - vertices(0,0)* vertices(2,1);
    //ak = vertices(0,0)* vertices(1,1) - vertices(1,0)* vertices(0,1);

    bi = vertices(1,1) - vertices(2,1);
    bj = vertices(2,1) - vertices(0,1);
    bk = vertices(0,1) - vertices(1,1);
    
    ci = vertices(2,0) - vertices(1,0);
    cj = vertices(0,0) - vertices(2,0);
    ck = vertices(1,0) - vertices(0,0);

    B(0,0) = 0.5 * bi / area;
    B(0,1) = 0.5 * bj / area;
    B(0,2) = 0.5 * bk / area;
    B(1,0) = 0.5 * ci / area;
    B(1,1) = 0.5 * cj / area;
    B(1,2) = 0.5 * ck / area;
}

void Tria3ThermalCalculator::ComputeStiffness
(
    const Matrix& vertices,    
    const ThermalMaterial& material,
    const Real& thickness,
    _OUT Matrix& ke
)
{
    Matrix D(s_nDim, s_nDim);
    Matrix B(s_nDim, s_nNode);
    Real area;
    Vector3D normal;
    
    const Real& k   = material.m_K;
    const Real& kx  = k;  //热传导系数
    const Real& ky  = k;  //热传导系数
    D(0, 0) = kx;
    D(1, 1) = ky;
    
    area = calc2DTriArea(vertices, normal);
    ComputeStrainMatrix(vertices, area, B);
    auto BT = B.Transpose ();
    ke  = BT * D * B * area * thickness;  //单元传热矩阵
}

void Tria3ThermalCalculator::ComputeCapacityStiffness
    (
    const SG::Algebra::Matrix& vertices, 
    const SG::DataStructure::FEM::ThermalMaterial& material,
    const Real& thickness,
    _OUT SG::Algebra::Matrix& ke3
    )
{
    Real area;
    Vector3D normal;
    area = calcTriArea (vertices, normal);
    
    ke3 = Matrix (3, 3, {2.0, 1.0, 1.0, 1.0, 2.0, 1.0, 1.0, 1.0, 2.0});
    const Real rho = material.m_RHO;
    const Real c = material.m_CP;
    ke3 = ke3 * (thickness * rho * c * area / 12.0);
}

void Tria3ThermalCalculator::ComputeHeatSourceLoad
    (
    const SG::Algebra::Matrix& coords,
    const std::vector<Real>& q,
    const Real& thickness,
    _OUT SG::Algebra::Matrix& p1
    )
{
    Real area;
    Vector3D normal;
    area = calcTriArea (coords, normal);
    Real qsum = (q[0] + q[1] + q[2]) / 3.0;

    p1 = Matrix (3, 1, {1.0, 1.0, 1.0});
    p1 = p1 * (thickness * qsum * area / 3.0);
}

void Tria3ThermalCalculator::ComputeElementMatries (const SG::DBManager::DBServiceFactorySP dbServiceFactory)
{
    auto& assembly = AssemblyEigen::GetInstance ();
    
    static const DofSet DOFs {DOFDirections::DOF_Ux};

    SG::DataStructure::FEM::CTRIA3Data eleData;
    PSHELLData property;
    ThermalMaterial material;
    
    std::vector<std::size_t> nodeIndexes(s_nNode);
    std::vector<Id> dispCoordIds(s_nNode);
    std::vector<int> dofGlobalIndex(s_nNode);
    std::vector<Point> points(s_nNode);
    
    const auto& tria3Service = dbServiceFactory ->get<IElementService>();
    const auto& nodeService = dbServiceFactory ->get<INodeService>();
    const auto eleSize = tria3Service->GetElementSize (m_type);
    for (std::size_t index{0}; index < eleSize; ++index)
    {
        /** 从数据库获取对应的单元数据 */
        tria3Service->FindByIndex (index, _OUT eleData);

        // Node coordinate
        ModelInfoQuery::QueryNode (nodeService, eleData.m_g, s_nNode, _OUT points, _OUT nodeIndexes, _OUT dispCoordIds);

        assembly.GetGlobalDOFId (nodeIndexes, DOFs, _OUT dofGlobalIndex);
        Matrix vertices(s_nNode,s_nDim);
        for (std::size_t iNode = 0; iNode < s_nNode; ++iNode)
        {
            vertices(iNode, 0) = points[iNode].m_x;
            vertices(iNode, 1) = points[iNode].m_y;
        }
         
        //获取材料参数
        dbServiceFactory ->get<IPSHELLService>()->FindById (eleData.m_pId, _OUT property);
        dbServiceFactory ->get<IThermalMaterialService>()->FindById (property.m_mId1, _OUT material);
        
        Matrix ke1(s_nNode, s_nNode);
        ComputeStiffness (vertices, material, property.m_thickness, _OUT ke1);
        assembly.AccumulateStiffness (dofGlobalIndex, ke1);
    }
}

// 后处理数据
void Tria3ThermalCalculator::DataRecovery (const SG::DBManager::DBServiceFactorySP dbServiceFactory, int subcaseId, SG::DataStructure::Common::Real time, const std::vector<SG::DataStructure::Common::Real>& globalDisplacement)
{
}

// 计算全局坐标系下单元应力刚度矩阵
void Tria3ThermalCalculator::ComputeStressStiffness (const SG::DBManager::DBServiceFactorySP dbServiceFactory, int staticSubId)
{
}
